# Replication archive for:
# Galos and Coppock. 2023. 
# Gender Composition Predicts Gender Bias: A Meta-reanalysis of Hiring Discrimination Audit Experiments. 
# Science Advances. Forthcoming.

library(tidyverse)
library(metafor)
library(broom)
source("code/helpers.R")

ates <- read_rds("data/estimated_ates.rds")
ates |> 
  ungroup() |> 
  summarize(n = n(),
            n_pos = sum(estimate >= 0))

field_meta_fit <-
  rma.uni(
    yi = estimate,
    sei = std.error,
    method = "REML",
    data = ates
  )


ates$weight <- weights(field_meta_fit)

field_meta_ests <- 
  tidy(field_meta_fit, conf.int = TRUE) %>%
  mutate(study_id = "Meta-analysis",
         study_name_citation = "Meta-analysis",
         facet = "Meta",
         weight = mean(ates$weight))


# Join study level with meta analysis -----------------------

gg_field <- 
  bind_rows(ates, field_meta_ests) %>% 
  ungroup() %>% 
  mutate(study_name_fct = fct_reorder(study_name_citation, estimate),
         entry = make_se_entry(estimate, std.error, 1),
         facet = factor(facet, levels = c("Study estimates", "Meta")))

# plot --------------------------------------------------------------------
g1 <- 
ggplot(gg_field, aes(estimate, study_name_fct)) +
  geom_point(aes(size = weight)) +
  geom_text(x = -22, aes(label = entry), hjust = 1, size = 2.5) +
  geom_linerange(aes(xmin = conf.low, xmax = conf.high), alpha = 0.3) +
  geom_vline(xintercept = 0, alpha = 0.25, linetype = "dashed") +
  geom_vline(data = select(field_meta_ests, -facet), aes(xintercept = estimate)) +
  geom_rect(
    data = select(field_meta_ests, -facet),
    aes(
      xmin = conf.low,
      xmax = conf.high,
      x = NULL,
      y = NULL
    ),
    ymin = 0,
    ymax = 100,
    alpha = .1
  ) +
  scale_size_area(max_size = 3) +
  theme_bw() +
  facet_grid(rows = vars(facet), scales = "free_y", space = "free") +
  coord_cartesian(xlim = c(-25, 25)) +
  theme(strip.background = element_blank(),
        strip.text = element_blank(),
        panel.grid = element_line(color = gray(0.95)),
        panel.grid.minor = element_blank(),
        axis.title.y = element_blank(),
        axis.text.y = element_text(size = 6),
        legend.position = "none") +
  labs(x = "Estimated average treatment effect in percentage points")

g1
